rm(list=ls())
# setwd(...) # set to appropriate working directory
data <- read.csv(file="MainData.csv")

# install necessary packages: install.packages("package name")
library(blockTools)
library(MatchIt)
library(aod)
library(sandwich)
library(stargazer)
library(lmtest)
library(MASS)
library(ordinal)
library(multiwayvcov)
library(car)
library(rms)
library(lme4)

#########################################################################################################################
### Table SI2.2.1: Descriptive Statistics ###
N <- as.data.frame(apply(data[,c(3:8, 11:13, 10)], 2, length)) # Number of observations
colnames(N) <- "N"
Mean <- as.data.frame(apply(data[,c(3:8, 11:13, 10)], 2, mean)) # Mean
colnames(Mean) <- "Mean"
SD <- as.data.frame(apply(data[,c(3:8, 11:13, 10)], 2, sd)) # St. dev.
colnames(SD) <- "St. dev."
Min <- as.data.frame(apply(data[,c(3:8, 11:13, 10)], 2, min))
colnames(Min) <- "Min."
Max <- as.data.frame(apply(data[,c(3:8, 11:13, 10)], 2, max))
colnames(Max) <- "Max."
cbind(N, Mean, SD, Min, Max)
#########################################################################################################################

#########################################################################################################################
### Figure SI2.2.1: Variable Densities ###
par(mfrow=c(3,3))
plot(density(data$Ethnic_Vote_Share), main='Ethnic Vote Share Density', cex.main=.8)
plot(density(data$Casualty), main='Casualty Density', cex.main=.8)
plot(density(data$Log_Casualty), main='Logged Casualty Density', cex.main=.8)
plot(density(data$HHI_SB), main='Serb-Bosniak HHI Density', cex.main=.8)
plot(density(data$HHI_CB), main='Croat-Bosniak HHI Density', cex.main=.8)
plot(density(data$HHI_SC), main='Serb-Croat HHI Density', cex.main=.8)
plot(density(data$Income_PC), main='Income Per Capita Density', cex.main=.8)
plot(density(data$Avg_Age), main='Average Age Density', cex.main=.8)
plot(density(data$Avg_Age_Sqr), main='Average Age Squared Density', cex.main=.8)
#########################################################################################################################

#########################################################################################################################
### Figure SI2.2.2: Variables Densities, Dichotomized ###
par(mfrow=c(3,3))
data_treated <- data[data$Casualty > mean(data$Casualty),]
data_untreated <- data[data$Casualty <= mean(data$Casualty),]

wilcox.test(data_treated$Ethnic_Vote_Share[data_treated$Year==1990], data_untreated$Ethnic_Vote_Share[data_untreated$Year==1990]) # not significant
# p = 0.108
plot(density(data_treated$Ethnic_Vote_Share[data_treated$Year==1990]), xlim=c(40,120), main='1990 Ethnic Vote Share Density', cex.main=.8, xlab='p = 0.108')
par(new=TRUE)
plot(density(data_untreated$Ethnic_Vote_Share[data_untreated$Year==1990]), xlim=c(40,120), lty=3, yaxt='n',xaxt='n', main='', xlab='')

wilcox.test(data_treated$HHI_SB, data_untreated$HHI_SB) 
# p = 0.074
plot(density(data_treated$HHI_SB), main='Serb-Bosniak HHI Density', cex.main=.8, xlab='p = 0.074', xlim=c(0,1))
par(new=TRUE)
plot(density(data_untreated$HHI_SB), main='', cex.main=.8, xlab='', xlim=c(0,1), lty=3, yaxt='n', xaxt='n')

wilcox.test(data_treated$HHI_CB, data_untreated$HHI_CB)
# p = 0.068
plot(density(data_treated$HHI_CB), main='Croat-Bosniak HHI Density', cex.main=.8, xlab='p = 0.068', xlim=c(0,1))
par(new=TRUE)
plot(density(data_untreated$HHI_CB), main='', cex.main=.8, xlab='', xlim=c(0,1), lty=3, yaxt='n', xaxt='n')

wilcox.test(data_treated$HHI_SC, data_untreated$HHI_SC) 
# p = 0.900
plot(density(data_treated$HHI_SC), main='Serb-Croat HHI Density', cex.main=.8, xlab='p = 0.900', xlim=c(0,1))
par(new=TRUE)
plot(density(data_untreated$HHI_SC), main='', cex.main=.8, xlab='', xlim=c(0,1), lty=3, yaxt='n', xaxt='n')

wilcox.test(data_treated$Income_PC, data_untreated$Income_PC)
# p = 0.011
plot(density(data_treated$Income_PC), main='Income Per Capita Density', cex.main=.8, xlab='p = 0.011', xlim=c(0,70000))
par(new=TRUE)
plot(density(data_untreated$Income_PC), main='', cex.main=.8, xlab='', xlim=c(0,70000), lty=3, yaxt='n', xaxt='n')

wilcox.test(data_treated$Avg_Age, data_untreated$Avg_Age)
# p = 0.086
plot(density(data_treated$Avg_Age), main='Average Age Density', cex.main=.8, xlab='p = 0.086', xlim=c(25,40))
par(new=TRUE)
plot(density(data_untreated$Avg_Age), main='', cex.main=.8, xlab='', xlim=c(25,40), lty=3, yaxt='n', xaxt='n')

wilcox.test(data_treated$Avg_Age_Sqr, data_untreated$Avg_Age_Sqr)
# p = 0.086
plot(density(data_treated$Avg_Age_Sqr), main='Average Age Squared Density', cex.main=.8, xlab='p = 0.086', xlim=c(600,1600))
par(new=TRUE)
plot(density(data_untreated$Avg_Age_Sqr), main='', cex.main=.8, xlab='', xlim=c(600,1600), lty=3, yaxt='n', xaxt='n')

plot(NULL)
legend('center', lty=c(1,3), legend = c('High Levels', 'Low Levels'), bty='n')
#########################################################################################################################

#########################################################################################################################
### Figure SI3.1: Model Diagnostic Plots ###
dummy.matrix <- as.matrix(cbind(as.numeric(data$Year == 2006), as.numeric(data$Year==2010), as.numeric(data$Year==2014)))

Log_Casualty.Model <- lm(Ethnic_Vote_Share ~ Log_Casualty:dummy.matrix+dummy.matrix+
                           Municipality-1, data)
par(mfrow=c(1,2), family='Times')
par(mfrow=c(1,2), family='Times')

# QQ-Plot
qqPlot(Log_Casualty.Model, col.lines='black', ylab = 'Studentized Residuals', main = 'QQ-Plot', cex.main=.4)

# Pearson Residuals on Fitted Values
plot(predict(Log_Casualty.Model), residuals(Log_Casualty.Model, type='pearson'), 
     xlab='Fitted Values', ylab='Pearson Residuals', 
     main='Pearson Residuals on Fitted Values',cex.main=1)
abline(h=0)
#########################################################################################################################

#########################################################################################################################
### Table SI4.1: Without Statistical Outliers ### 
# Model 1
outlierTest(Log_Casualty.Model, cutoff=0.05)
data[c(322, 79, 262), c("Municipality", "Year")]
outly_data <- data[-c(79,262),]
outly_dummy <- dummy.matrix[-c(79,262),]
Log_Casualty.Outliers <- lm(Ethnic_Vote_Share ~ Log_Casualty:outly_dummy
                            + outly_dummy+Municipality-1, outly_data)

vcovCL <- cluster.vcov(Log_Casualty.Outliers, ~outly_data$Municipality+outly_data$Year)
coeftest(Log_Casualty.Outliers, vcovCL) # results
nobs(Log_Casualty.Outliers) # 426 observations
summary(Log_Casualty.Outliers)$r.squared # 0.9925671

# Model 2 
which(data$Municipality == "Velika Kladusa")
outly_data2 <- data[-c(79,262,1,2,3,322),]
outly_dummy2 <- dummy.matrix[-c(79,262,1,2,3,322),]
Log_Casualty.Outliers2 <- lm(Ethnic_Vote_Share ~ Log_Casualty:outly_dummy2
                            + outly_dummy2+Municipality-1, outly_data2)

vcovCL <- cluster.vcov(Log_Casualty.Outliers2, ~outly_data2$Municipality+outly_data2$Year)
coeftest(Log_Casualty.Outliers2, vcovCL) # results
nobs(Log_Casualty.Outliers2) # 422 observations
summary(Log_Casualty.Outliers2)$r.squared # 0.9938933
#########################################################################################################################

#########################################################################################################################
### Table SI4.2: Without Ethnic Enclaves ###
enclave_data <- data[-c(79,82,262),]
enclave_dummy <- dummy.matrix[-c(79,82,262),]
Log_Casualty.Enclaves <- lm(Ethnic_Vote_Share ~ Log_Casualty:enclave_dummy
                         + enclave_dummy+Municipality-1, enclave_data)

vcovCL <- cluster.vcov(Log_Casualty.Enclaves, ~enclave_data$Municipality+enclave_data$Year)
coeftest(Log_Casualty.Enclaves, vcovCL) # results
nobs(Log_Casualty.Enclaves) # 425 observations
summary(Log_Casualty.Enclaves)$r.squared # 0.9927555
#########################################################################################################################

#########################################################################################################################
### Table SI4.3: Casualty Unlogged ###
Casualty.Model <- lm(Ethnic_Vote_Share ~ Casualty:dummy.matrix
                            + dummy.matrix+Municipality-1, data)
vcovCL <- cluster.vcov(Casualty.Model, ~data$Municipality+data$Year)
coeftest(Casualty.Model, vcovCL) # results
nobs(Casualty.Model) # 428 observations
summary(Casualty.Model)$r.squared # 0.9905516
#########################################################################################################################

#########################################################################################################################
### Table SI4.4: Civlian Casualties ###
Civilian.Model <- lm(Ethnic_Vote_Share ~ CivilianCasualty:dummy.matrix+dummy.matrix+
                       Municipality-1, data)
vcovCL <- cluster.vcov(Civilian.Model, ~data$Municipality+data$Year)
coeftest(Civilian.Model, vcovCL) # results
nobs(Civilian.Model) # 404 observations
summary(Civilian.Model)$r.squared # 0.9904033
#########################################################################################################################

#########################################################################################################################
### Table SI4.5: Null Model ###
Null.Model <- lm(Ethnic_Vote_Share ~ dummy.matrix+
                        Municipality-1, data)

vcovCL <- cluster.vcov(Null.Model, ~data$Municipality+data$Year)
coeftest(Null.Model, vcovCL) # results
nobs(Null.Model) # 428 observations
summary(Null.Model)$r.squared # 0.990412

vcovCL <- cluster.vcov(Log_Casualty.Model, ~data$Municipality+data$Year)

# Wald Test comparing null model fit with main model fit. The main model performs significantly better than the null model.
waldtest(Null.Model, Log_Casualty.Model, vcov=vcovCL) 
#########################################################################################################################

#########################################################################################################################
### Table SI4.6: Dichotomized Treatment ###
treated <- data$Casualty>mean(data$Casualty)
Log_Casaulty.Dich <- lm(Ethnic_Vote_Share ~ treated:dummy.matrix+dummy.matrix+
                        Municipality-1, data)
vcovCL <- cluster.vcov(Log_Casaulty.Dich, ~data$Municipality+data$Year)
coeftest(Log_Casaulty.Dich, vcovCL) # results
nobs(Log_Casaulty.Dich) # 428 observations
summary(Log_Casaulty.Dich)$r.squared # 0.9906337
#########################################################################################################################

#########################################################################################################################
### Table SI4.7: Unpartitioned Municipalities ###
index <- which(data$Partition == 0)
UnpartitionedData <- data[index,]
UnpartitionedDummy <- dummy.matrix[index,]

Unpartitioned.Model <- lm(Ethnic_Vote_Share ~ Log_Casualty:UnpartitionedDummy + UnpartitionedDummy +
                            Municipality-1, UnpartitionedData)
vcovCL <- cluster.vcov(Unpartitioned.Model, ~UnpartitionedData$Municipality+UnpartitionedData$Year)
coeftest(Unpartitioned.Model, vcovCL) # results
nobs(Unpartitioned.Model) # 304 observations
summary(Unpartitioned.Model)$r.squared # 0.9907471
#########################################################################################################################

#########################################################################################################################
### Table SI4.8: Mediating Effects of Homogeneity ###
rm(list=ls())
data <- read.csv(file="MainData.csv")
HHI.2013_Data <- data[complete.cases(data$HHI.2013),]
matrix.2013Y <- as.matrix(cbind(as.numeric(HHI.2013_Data$Year == 2006), as.numeric(HHI.2013_Data$Year==2010), as.numeric(HHI.2013_Data$Year==2014)))

first <- lmer(Ethnic_Vote_Share ~ Log_Casualty:matrix.2013Y[,1] 
              + Log_Casualty:matrix.2013Y[,2] 
              + Log_Casualty:matrix.2013Y[,3]
              + HHI_SB
              + HHI_CB
              + HHI_SC
              + HHI.2013
              + Pop_Density
              + Income_PC
              + Avg_Age
              + Avg_Age_Sqr 
              + matrix.2013Y + (1|Municipality), HHI.2013_Data)

direct <- lmer(I(Ethnic_Vote_Share - fixef(first)['HHI.2013']*HHI.2013) ~ Log_Casualty:matrix.2013Y[,1]
               + Log_Casualty:matrix.2013Y[,2] 
               + Log_Casualty:matrix.2013Y[,3]
               + HHI_SB
               + HHI_CB
               + HHI_SC
               + Pop_Density
               + Income_PC
               + Avg_Age
               + Avg_Age_Sqr 
               + matrix.2013Y + (1|Municipality), HHI.2013_Data)

boots <- 1000
fl.boots <- matrix(NA, nrow=boots, ncol=14)
for(b in 1:boots){
  d.star <- HHI.2013_Data[sample(1:nrow(HHI.2013_Data), replace=TRUE),]
  matrix.2013Y <- as.matrix(cbind(as.numeric(d.star$Year == 2006), as.numeric(d.star$Year==2010), as.numeric(d.star$Year==2014)))
  first.boots <- lmer(Ethnic_Vote_Share ~ Log_Casualty:matrix.2013Y[,1] 
                      + Log_Casualty:matrix.2013Y[,2] 
                      + Log_Casualty:matrix.2013Y[,3]
                      + HHI_SB
                      + HHI_CB
                      + HHI_SC
                      + HHI.2013
                      + Pop_Density
                      + Income_PC
                      + Avg_Age
                      + Avg_Age_Sqr 
                      + matrix.2013Y + (1|Municipality), d.star)
  
  direct.boots <- lmer(I(Ethnic_Vote_Share - fixef(first)['HHI.2013']*HHI.2013) ~ Log_Casualty:matrix.2013Y[,1]
                       + Log_Casualty:matrix.2013Y[,2] 
                       + Log_Casualty:matrix.2013Y[,3]
                       + HHI_SB
                       + HHI_CB
                       + HHI_SC
                       + Pop_Density
                       + Income_PC
                       + Avg_Age
                       + Avg_Age_Sqr 
                       + matrix.2013Y + (1|Municipality), d.star)
  fl.boots[b,] <- fixef(direct.boots)
  
}
fixef(direct) # coefficients

# The standard errors are bootstrapped, and will therefore vary slightly with each iteration.
(SEs <- apply(fl.boots,2,sd)) # standard erros

# Assessing significance with 95% confidence intervals
rep(fixef(direct),each=2)+c(-1,1)*1.96*rep(SEs,each=2)
# 2006, 2010, and 2014 significant

nobs(direct) # 412 observations
AIC(direct) # 3155.452
#########################################################################################################################

#########################################################################################################################
### Table SI4.9: Mediating Effects of Homogeneity Over Time ###
rm(list=ls())
data <- read.csv(file="MainData.csv")
TimedHHI_Data <- data[complete.cases(data$TimedHHI),]
Timed.dummy <- as.matrix(cbind(as.numeric(TimedHHI_Data$Year == 2006), as.numeric(TimedHHI_Data$Year==2010), as.numeric(TimedHHI_Data$Year==2014)))

first <- lm(Ethnic_Vote_Share ~ Log_Casualty:Timed.dummy[,1] 
            + Log_Casualty:Timed.dummy[,2] 
            + Log_Casualty:Timed.dummy[,3]
            + Timed.dummy + Municipality-1+TimedHHI, TimedHHI_Data)

direct <- lm(I(Ethnic_Vote_Share - coef(first)['TimedHHI']*TimedHHI) ~ Log_Casualty:Timed.dummy[,1]
             + Log_Casualty:Timed.dummy[,2] 
             + Log_Casualty:Timed.dummy[,3]
             + Timed.dummy + Municipality-1, TimedHHI_Data)

boots <- 1000
fl.boots <- matrix(NA, nrow=boots, ncol=6)
for(b in 1:boots){
  d.star <- TimedHHI_Data[sample(1:nrow(TimedHHI_Data), replace=TRUE),]
  Timed.dummy <- as.matrix(cbind(as.numeric(d.star$Year == 2006), as.numeric(d.star$Year==2010), as.numeric(d.star$Year==2014)))
  first.boots <- lm(Ethnic_Vote_Share ~ Log_Casualty:Timed.dummy[,1] 
                    + Log_Casualty:Timed.dummy[,2] 
                    + Log_Casualty:Timed.dummy[,3]
                    + TimedHHI
                    + Timed.dummy + Municipality-1, d.star)
  
  direct.boots <- lm(I(Ethnic_Vote_Share - coef(first)['TimedHHI']*TimedHHI) ~ Log_Casualty:Timed.dummy[,1]
                     + Log_Casualty:Timed.dummy[,2] 
                     + Log_Casualty:Timed.dummy[,3]
                     + Timed.dummy + Municipality-1, d.star)
  
  fl.boots[b,] <- coef(direct.boots)[c(1:3, (length(coef(direct.boots))-2):length(coef(direct.boots)))]
  
}
coef(direct)[c(1:3, 91:93)] # coefficients

# The standard errors are bootstrapped, and will therefore vary slightly with each iteration.
(SEs <- apply(fl.boots,2,sd)) # standard errors

# Assessing significance with 95% confidence intervals
rep(coef(direct)[c(1:3, (length(coef(direct))-2):length(coef(direct)))],each=2)+c(-1,1)*1.96*rep(SEs,each=2) # 95% (p < 0.05)
# 2006 and 2010 significant

nobs(direct) # 348 observations
summary(direct)$r.squared # 0.9883608
#########################################################################################################################

#########################################################################################################################
### Table SI4.10: Mediating Effects of Population Displacement ###
# Model 1: Net Displacement
rm(list=ls())
data <- read.csv(file="MainData.csv")
ND_Data <- data[complete.cases(data$Net_Displacement),]
ND.matrix <- as.matrix(cbind(as.numeric(ND_Data$Year == 2006), as.numeric(ND_Data$Year==2010), as.numeric(ND_Data$Year==2014)))

first <- lm(Ethnic_Vote_Share ~ Log_Casualty:ND.matrix[,1] 
            + Log_Casualty:ND.matrix[,2] 
            + Log_Casualty:ND.matrix[,3]
            + ND.matrix + Municipality-1+Net_Displacement, ND_Data)

direct <- lm(I(Ethnic_Vote_Share - coef(first)['Net_Displacement']*Net_Displacement) ~ Log_Casualty:ND.matrix[,1]
             + Log_Casualty:ND.matrix[,2] 
             + Log_Casualty:ND.matrix[,3]
             + ND.matrix + Municipality-1, ND_Data)

boots <- 1000
fl.boots <- matrix(NA, nrow=boots, ncol=6)
for(b in 1:boots){
  d.star <- ND_Data[sample(1:nrow(ND_Data), replace=TRUE),]
  ND.matrix <- as.matrix(cbind(as.numeric(d.star$Year == 2006), as.numeric(d.star$Year==2010), as.numeric(d.star$Year==2014)))
  first.boots <- lm(Ethnic_Vote_Share ~ Log_Casualty:ND.matrix[,1] 
                    + Log_Casualty:ND.matrix[,2] 
                    + Log_Casualty:ND.matrix[,3]
                    + Net_Displacement
                    + ND.matrix + Municipality-1, d.star)
  
  direct.boots <- lm(I(Ethnic_Vote_Share - coef(first)['Net_Displacement']*Net_Displacement) ~ Log_Casualty:ND.matrix[,1]
                     + Log_Casualty:ND.matrix[,2] 
                     + Log_Casualty:ND.matrix[,3]
                     + ND.matrix + Municipality-1, d.star)
  
  fl.boots[b,] <- coef(direct.boots)[c(1:3, (length(coef(direct.boots))-2):length(coef(direct.boots)))]
  
}
coef(direct)[c(1:3, 107:109)] # coefficients

# The standard errors are bootstrapped, and will therefore vary slightly with each iteration.
(SEs <- apply(fl.boots,2,sd)) # standard errors

# Assessing significance with 95% confidence intervals
rep(coef(direct)[c(1:3, (length(coef(direct))-2):length(coef(direct)))],each=2)+c(-1,1)*1.96*rep(SEs,each=2) # 95% (p < 0.05)
# 2006 and 2010 significant 

nobs(direct) # 412 observations
summary(direct)$r.squared # 0.9909547

# Model 2: Gross Displacement
rm(list=ls())
data <- read.csv(file="MainData.csv")
GD_Data <- data[complete.cases(data$Gross_Displacement),]
GD.matrix <- as.matrix(cbind(as.numeric(GD_Data$Year == 2006), as.numeric(GD_Data$Year==2010), as.numeric(GD_Data$Year==2014)))

first <- lm(Ethnic_Vote_Share ~ Log_Casualty:GD.matrix[,1] 
            + Log_Casualty:GD.matrix[,2] 
            + Log_Casualty:GD.matrix[,3]
            + GD.matrix + Municipality-1+Gross_Displacement, GD_Data)

direct <- lm(I(Ethnic_Vote_Share - coef(first)['Gross_Displacement']*Gross_Displacement) ~ Log_Casualty:GD.matrix[,1]
             + Log_Casualty:GD.matrix[,2] 
             + Log_Casualty:GD.matrix[,3]
             + GD.matrix + Municipality-1, GD_Data)

boots <- 1000
fl.boots <- matrix(NA, nrow=boots, ncol=6)
for(b in 1:boots){
  d.star <- GD_Data[sample(1:nrow(GD_Data), replace=TRUE),]
  GD.matrix <- as.matrix(cbind(as.numeric(d.star$Year == 2006), as.numeric(d.star$Year==2010), as.numeric(d.star$Year==2014)))
  first.boots <- lm(Ethnic_Vote_Share ~ Log_Casualty:GD.matrix[,1] 
                    + Log_Casualty:GD.matrix[,2] 
                    + Log_Casualty:GD.matrix[,3]
                    + Gross_Displacement
                    + GD.matrix + Municipality-1, d.star)
  
  direct.boots <- lm(I(Ethnic_Vote_Share - coef(first)['Gross_Displacement']*Gross_Displacement) ~ Log_Casualty:GD.matrix[,1]
                     + Log_Casualty:GD.matrix[,2] 
                     + Log_Casualty:GD.matrix[,3]
                     + GD.matrix + Municipality-1, d.star)
  
  fl.boots[b,] <- coef(direct.boots)[c(1:3, (length(coef(direct.boots))-2):length(coef(direct.boots)))]
  
}
coef(direct)[c(1:3, 107:109)] # coefficients

# The standard errors are bootstrapped, and will therefore vary slightly with each iteration.
(SEs <- apply(fl.boots,2,sd)) # standard errors

# Assessing significance with 95% confidence intervals
rep(coef(direct)[c(1:3, (length(coef(direct))-2):length(coef(direct)))],each=2)+c(-1,1)*1.96*rep(SEs,each=2) # 95% (p < 0.05)
# 2006 and 2010 significant

nobs(direct) # 412 observations
summary(direct)$r.squared # 0.9911882
#########################################################################################################################

#########################################################################################################################
### Table SI4.11 ###
# Model 1: log(Refugees)
rm(list=ls())
data <- read.csv(file="MainData.csv")

dummy.matrix <- as.matrix(cbind(as.numeric(data$Year == 2006), as.numeric(data$Year==2010), as.numeric(data$Year==2014)))

Log_Refugees.Model <- lm(Ethnic_Vote_Share ~ log(Refugees):dummy.matrix + dummy.matrix +
                            Municipality-1, data)
vcovCL <- cluster.vcov(Log_Refugees.Model, ~data$Municipality+data$Year)
coeftest(Log_Refugees.Model, vcovCL) # results
nobs(Log_Refugees.Model) # 428 observations
summary(Log_Refugees.Model)$r.squared # 0.9910958

# Model 2: log(Prison Camps) 
Log_Prison_Camps.Model <- lm(Ethnic_Vote_Share ~ log(Prison_Camps+0.001):dummy.matrix + dummy.matrix +
                           Municipality-1, data)
vcovCL <- cluster.vcov(Log_Prison_Camps.Model, ~data$Municipality+data$Year)
coeftest(Log_Prison_Camps.Model, vcovCL) # results
nobs(Log_Prison_Camps.Model) # 388 observations
summary(Log_Prison_Camps.Model)$r.squared # 0.9908079
#########################################################################################################################

#########################################################################################################################
# Survey Data
SurveyData <- read.csv("SurveyData.csv")

Ethnic_Friends_Data <- na.omit(SurveyData[,c('Ethnic_Friends','Violence','Gender','Age','HHI_SB','HHI_CB','HHI_SC',
                                            'Pop_Density','Income_PC','Municipality')])
Closest_Friends_Data <- na.omit(SurveyData[,c('Closest_Friends','Violence','Gender','Age','HHI_SB','HHI_CB','HHI_SC',
                                              'Pop_Density','Income_PC','Municipality')])
National_Trust_Data <- na.omit(SurveyData[,c('National_Trust','Violence','Gender','Age','HHI_SB','HHI_CB','HHI_SC',
                                             'Pop_Density','Income_PC','Municipality')])
Representation_Data <- na.omit(SurveyData[,c('Representation','Violence','Gender','Age','HHI_SB','HHI_CB','HHI_SC',
                                             'Pop_Density','Income_PC','Municipality')])

# For Table SI4.12, we present the code for the right columns of each model. These are the models
# where we use the self-reported measure (Violence).

# The left column of each model are the results from Table 2 in the main text. To see the code for those
# models, please see replication material for content presented in the main text.

### Table SI4.12 ###
# Model 1: Ethnic Friends, Right Column
Ethnic_Friends.model <- lrm(Ethnic_Friends ~ Violence
                            + Gender
                            + Age
                            + HHI_SB
                            + HHI_CB
                            + HHI_SC
                            + Pop_Density
                            + Income_PC, y=TRUE, x=TRUE, data=Ethnic_Friends_Data)
robcov(Ethnic_Friends.model,cluster=Ethnic_Friends_Data$Municipality)

# Model 2: Closest Friends, Right Column
Closest_Friends.model <- lrm(Closest_Friends ~ Violence
                             + Gender
                             + Age
                             + HHI_SB
                             + HHI_CB
                             + HHI_SC
                             + Pop_Density
                             + Income_PC, y=TRUE, x=TRUE, data=Closest_Friends_Data)
robcov(Closest_Friends.model,cluster=Closest_Friends_Data$Municipality)

# Model 3: National Trust, Right Column ###
National_Trust.model <- lrm(National_Trust ~ Violence
                            + Gender
                            + Age
                            + HHI_SB
                            + HHI_CB
                            + HHI_SC
                            + Pop_Density
                            + Income_PC, y=TRUE, x=TRUE, data=National_Trust_Data)
robcov(National_Trust.model,cluster=National_Trust_Data$Municipality)

# Model 4: Representation, Right Column ###
Representation.model <- lrm(Representation ~ Violence
                            + Gender
                            + Age
                            + HHI_SB
                            + HHI_CB
                            + HHI_SC
                            + Pop_Density
                            + Income_PC, y=TRUE, x=TRUE, data=Representation_Data)
robcov(Representation.model,cluster=Representation_Data$Municipality)
#########################################################################################################################

#########################################################################################################################
### Table SI5.1: Nearest Neighbor Matching ###
matchData <- read.csv(file="Matching.csv")
median(matchData$Casualty) # 1.766209

for(i in 1:nrow(matchData)){
  if(matchData$Casualty[i] < 1.766209){
    matchData$Treat[i] <- 0
  }
  if(matchData$Casualty[i] >= 1.766209){
    matchData$Treat[i] <- 1
  }
}

match.out <- matchit(Treat ~ Pop_Density + Income_PC + HHI_BS + HHI_BC + HHI_SC + Avg_Age + Avg_Age_Sqr,
                     data=matchData, method="nearest", distance="mahalanobis", replace=TRUE)
(matched.cases <- match.out$match.matrix) # First column treated, second column control.

# Vitez and Kakanj is the matched pair we select for comparison. This is pair 44, 33.
matchData[44,] # Vitez
matchData[33,] # Kakanj

K <- data[data$Municipality == "Kakanj",]
V <- data[data$Municipality == "Vitez",]

K <- c(unique(K$Casualty), unique(K$Pop_Density), unique(K$Income_PC), unique(K$Avg_Age), unique(K$Avg_Age_Sqr),
       unique(K$HHI_SB), unique(K$HHI_CB), unique(K$HHI_SC), 
       K$Ethnic_Vote_Share[4], K$Ethnic_Vote_Share[1], K$Ethnic_Vote_Share[2], K$Ethnic_Vote_Share[3],
       K$Ethnic_Vote_Share[1]-K$Ethnic_Vote_Share[4], K$Ethnic_Vote_Share[2]-K$Ethnic_Vote_Share[4], K$Ethnic_Vote_Share[3]-K$Ethnic_Vote_Share[4])

V <- c(unique(V$Casualty), unique(V$Pop_Density), unique(V$Income_PC), unique(V$Avg_Age), unique(V$Avg_Age_Sqr),
       unique(V$HHI_SB), unique(V$HHI_CB), unique(V$HHI_SC), 
       V$Ethnic_Vote_Share[4], V$Ethnic_Vote_Share[1], V$Ethnic_Vote_Share[2], V$Ethnic_Vote_Share[3],
       V$Ethnic_Vote_Share[1]-V$Ethnic_Vote_Share[4], V$Ethnic_Vote_Share[2]-V$Ethnic_Vote_Share[4], V$Ethnic_Vote_Share[3]-V$Ethnic_Vote_Share[4])

TSI5.1 <- cbind(K, V)
colnames(TSI5.1) <- c("Kakanj", "Vitez")
rownames(TSI5.1) <- c("Casualty", "Population density", "Income per capita", "Average age", "Average age squared",
                      "HHI Serb-Bosniak", "HHI Croat-Bosniak", "HHI Serb-Croat", "Ethnic Voting 1990", "Ethnic Voting 2006",
                      "Ethnic Voting 2010", "Ethnic Voting 2014", "Ethnic difference 2006", "Ethnic difference 2010",
                      "Ethnic difference 2014")
TSI5.1

#########################################################################################################################
### Table SI5.2: High Dimensional Blocking ###
blockData <- read.csv(file="Matching.csv")
(out <- block(blockData,n.tr = 2, 
              id.vars = c("Municipality"), 
              block.vars=c("Pop_Density", "Income_PC", "HHI_BS", "HHI_BC", "HHI_SC", "Avg_Age", "Avg_Age_Sqr"), 
              algorithm="optGreedy", 
              distance="mahalanobis"))
# Sanski Most and Prijedor selected for comparison. This pair is located in row 3.

S <- data[data$Municipality == "Sanski Most",]
P <- data[data$Municipality == "Prijedor",]

S <- c(unique(S$Casualty), unique(S$Pop_Density), unique(S$Income_PC), unique(S$Avg_Age), unique(S$Avg_Age_Sqr),
       unique(S$HHI_SB), unique(S$HHI_CB), unique(S$HHI_SC), 
       S$Ethnic_Vote_Share[4], S$Ethnic_Vote_Share[1], S$Ethnic_Vote_Share[2], S$Ethnic_Vote_Share[3],
       S$Ethnic_Vote_Share[1]-S$Ethnic_Vote_Share[4], S$Ethnic_Vote_Share[2]-S$Ethnic_Vote_Share[4], S$Ethnic_Vote_Share[3]-S$Ethnic_Vote_Share[4])

P <- c(unique(P$Casualty), unique(P$Pop_Density), unique(P$Income_PC), unique(P$Avg_Age), unique(P$Avg_Age_Sqr),
       unique(P$HHI_SB), unique(P$HHI_CB), unique(P$HHI_SC), 
       P$Ethnic_Vote_Share[4], P$Ethnic_Vote_Share[1], P$Ethnic_Vote_Share[2], P$Ethnic_Vote_Share[3],
       P$Ethnic_Vote_Share[1]-P$Ethnic_Vote_Share[4], P$Ethnic_Vote_Share[2]-P$Ethnic_Vote_Share[4], P$Ethnic_Vote_Share[3]-P$Ethnic_Vote_Share[4])

TSI5.2 <- cbind(S, P)
colnames(TSI5.2) <- c("Sanski Most", "Prijedor")
rownames(TSI5.2) <- c("Casualty", "Population density", "Income per capita", "Average age", "Average age squared",
                      "HHI Serb-Bosniak", "HHI Croat-Bosniak", "HHI Serb-Croat", "Ethnic Voting 1990", "Ethnic Voting 2006",
                      "Ethnic Voting 2010", "Ethnic Voting 2014", "Ethnic difference 2006", "Ethnic difference 2010",
                      "Ethnic difference 2014")
TSI5.2
#########################################################################################################################